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Semi-Spectral Method for the Wigner equation 
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We propose a numerical method to solve the Wigner equation in quantum systems of spinless, 
non-relativistic particles. The method uses a spectral decomposition into basis functions 

in momentum-space to obtain a system of first-order advection-reaction equations. The resulting 
equations are solved by splitting the reaction and advection steps so as to allow the combination 
of numerical techniques from quantum mechanics and computational fluid dynamics by identifying 
the skew-hermitian reaction matrix as a generator of unitary rotations. The method is validated for 
the case of particles subject to a one-dimensional (an-)harmonic potential using finite-differences for 
the advection part. Thereby, we verify the second order of convergence and observe non-classical 
behavior in the evolution of the Wigner function. 
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I. INTRODUCTION 

The Wigner formalisrm also known as quantum me¬ 
chanics in phase space [l|, provides an alternative but 
equivalent description of quantum mechanics in trems 
of a (quasi)-distribution function of the particle position 
and momentum. It has proven to provide a helpful sup¬ 
plement to operator methods in Hilbert space as well as 
to path integral formulations, and has offered new in¬ 
sights into the relation between quantum and classical 
physics, as it does not discriminate between coordinate 
and momentum space. For instance, it has been a fruitful 
perspective for the study of quantum chaos. In addition, 
it offers the opportunity to systematically consider quan¬ 
tum corrections to the classical dynamics by expanding 
the quantum Liouville equation around h pz 0. Nowadays, 
it is also a valuable tool in the fields of quantum optics 
as well as nuclear, plasma and semiconductor physics to 
describe transport processes, for example, in open quan¬ 
tum systems [l|. The Wigner function, introduced by 
E. Wigner in Ref. 0, is the Weyl transformation of the 
density matrix and a quasi-probability distribution that 
can “intuitively” account for scattering and decoherence 
effects in quantum systems Si- It differs from a classi¬ 
cal probability distributions as it can change sign during 
the evolution especially in regions where quantum inter¬ 
ference effects become important. 

Since the Wigner equation was introduced in 1932, it 
has been tackled by various numerical approaches, such 
as finite differences S i, Fourier spectral collocation 
3, B 1 deterministic particle S [13 , and Monte-Carlo 
III ll j . Here, we extend the technique described in 
Refsrp,S| to arbitrary basis functions (finip) of in 

momentum-space and reveal the underlying mathemati¬ 
cal structure of the resulting infinite-dimensional set of 
reaction-advection equations. By using this formulation 
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we show that the action of the potential on the Wigner 
function is a unitary rotation of its coefficient vector, 
whereas the advection operation can be discretized by 
various techniques used in computational fluid dynamics, 
such as finite difference, finite volume or finite element, 
cf. Ref. [l3- In that way, one is able to construct a 
finite element simulation of the Wigner evolution. Em¬ 
ploying a more general basis we assume that the higher 
computational costs of our method, compared to 

0{N\ogN) for the spectral Fourier decomposition are 
outweighed by a smaller number N of basis functions to 
obtain the same order of accuracy through focusing the 
computational effort to regions of interest, as for example 
in the case of Wigner functions that are strongly local¬ 
ized in momentum-space, such as particles in a periodic 
potential, cf. Bloch’s theorem and Ref. [l^ • In addition, 
the “artificial” periodization of the Wigner function can 
be avoided, which may mitigate the self-interaction of the 
distribution at the domain boundaries of the simulation 
that is present for the Fourier basis choice, cf. Ref. (T^ - 


This article is organized as follows. First, we give 
an introduction to the Wigner formalism and present 
the properties of the Wigner equation, especially for 
the pseudo-differential operator. In section IHII we show 
the details of the numerical method to handle the ob¬ 
tained multi-dimensional reaction-advection equation. 
Thirdly, we validate the technique by simulating a one¬ 
dimensional (an-)harmonic oscillator, which offers the 
opportunity to compare with an analytical solution, cf. 
Ref. [ 13 , such that we can perform a convergence analy¬ 
sis, and to observe quantum effects when the anharmonic 
potential is used. To study tunneling phenomena, we 
show the evolution of bounded states in the double well 
potential and measure the spread as well as the covari¬ 
ance of the Wigner function in phase space. In the last 
section, we will highlight the strengths and weaknesses 
of the approach. 
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II. WIGNER FORMALISM 


and the canonical commutation relation is written as 


Our aim is to simulate the time-evolution of the 
Wigner function w{T,q,p) of a d-dimensional system of 
non-relativistic spinless particles of mass m subject to 
the potential U{T,q), based on the Wigner equation 


0 = drW H-- V„w -I- Q[U]w , 

m 


( 1 ) 


0[t/] = - p(T, q + ihVp/2) - U{t, q- *Wp/2)J . (2) 


The independent variables are time r, space q and mo¬ 
mentum p respectively. Hence, the Wigner function itself 
has the dimension h~‘^, since it fulfills 


[xj,Vk] = leSj^k ■ ( 10 ) 

The symbol 0[H] stands for the pseudo-differential oper¬ 
ator, whose action on W can be written as an integral 

iB C 

Q[V]W = — df] 5V{t, X, rJ)W{t, x, , (11) 

R'i 

SV{t, x,fj) = V (t, x+^ff^ -V (t,x- I??) , (12) 

W{t,x,f]) = jdp W{t,x,p)e^^'P . (13) 


Jdq Jdp w{t, q,p)= Np , (3) 

R<i R'i 


where Np is the number of particles in the system. To 
have an easier grasp on the equation, we make it dimen¬ 
sionless. First of all, we measure the Wigner function 
with respect to i.e. we introduce the dimension¬ 
less Wigner function W{T,q,p) = hfw{T,q,p)/Np, and 
we employ the following scaling relations 

x = q/l ,t = TjT , V = ^p ,V{t,x) = U (t, q)/U , 

described in Ref. Thus, we obtain the dimensionless 
Wigner equation 


0 = dtW + v- -f <d[V]W , 


( 4 ) 


0[H] = f ^ 




where we have introduced the dimensionless constants 


_ hT UT^ 

Pm ’ Pm ’ 


( 6 ) 


which we call effective Planck’s constant and potential 
strength, respectively. The names are not arbitrary and 
reflect the natural occurrence of these numbers in the 
dimensionless formulation of the dynamics. For instance, 
Eq. ([U directly translates into 


^—d 




( 7 ) 


where we have replaced Planck’s constant by its scaled 
counterpart; the Wigner transform of a pure state 'I' be¬ 
comes 

W = Jiy *• (i. X + f) f (f, df - f) 9? , (8) 

R‘^ 


If the potential is locally well-approximated by a Taylor 
series around e ~ 0 we can write 

V{t,x + erf/2) « {e/2f^ (^ 4 ) 

|A|=0 

where A is a multi-index of dimension d, i.e. |A| = 

A! = and = ULidp There¬ 

fore, the action of the pseudo-differential operator on the 
Wigner function reads 

e[V]W = -B Yy {D^W) . 

|A|GNodd 

(15) 

Nodd stands for strictly positive odd integers, such that 
the sum is always real. For this treatment, the potential 
needs to be an analytic function defined on an open set 
D C X R, i.e. explicit dependence on time is possible. 
This implies two important properties for us. Firstly, the 
function is locally given by a convergent power or Taylor 
series. Secondly, one can find an upper bound for all 
derivatives of the function, since for every compact set 
K C D, for all (t,x) S K, and for all |A| G No there 
exists a constant C such that 

\D^V\ < Cl^l+^A! . ( 16 ) 

Hence, by choosing the right time and length scale (T, 1) 
for the problem one can find a convergent power series 
representation of the pseudo-differential operator. The 
correct choice means 

//^p\ 1^1 

(-) =“■ 

such that the series in Eq. (fT^ converges (locally) uni¬ 
formly under the assumption that D^W is bounded. 


where 'i’{t,x) = l'^^'^'ijj{T,q) is the dimensionless wave 
function in position space; the dimensionless time- 
dependent Schrodinger equation reads 




\d\'^/2 + BV{t,x) 


4^ , 


(9) 


A. General properties 

Expanding the Wigner function in momentum-space 
into a set of orthonormal basis functions of 
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with the inner product 


{(l>iA 3)2 = fdv (l)*{v)^j{v) = , (18) 

R‘^ 

meaning that 

= ^ak{t,x) (l)k{v) , (19) 

fcGN 

we can rewrite the Wigner equation, Eq. (SD, into an infi¬ 
nite system of linear, first-order partial differential equa¬ 
tions (PDEs) for the coefficients ak{t, x) G C. The system 
can be derived by using the orthonormality property of 
the basis functions, Eq. dni). Depending on the choice 
of the basis we will find different sets of PDEs. In gen¬ 
eral, all the sets can be written as a multi-dimensional 
reaction-advection equation 


d 

did -b A^’^^dxjd + Mv{t, x)d = 0 , ( 20 ) 

i=l 


where Ad\Mv{t,x) are square matrices and a = 
( 01 , 02 , 03 ,...) is the coefficient vector. Independent of 
the basis choice, the matrix Mv(t,x) is skew-hermitian, 
which can be demonstrated employing Eq. (ESI or Eq. 
(ED. When using formula ((T5|) we have to assume that 
the basis functions are Under this condition we 

can shift the uneven derivatives, |A| S Nodd, which ap¬ 
pear as summands in the pseudo-differential operator, to 
show the skew-hermiticity. Demonstrating this property 
for a general set of basis functions of L^(R‘^), i.e. even 
non-differentiable, we use Eq. ED to write 


(Mvd), = jdv ct^liv) (0[U]W) Xx,v) (21) 

R'i 

iB f 

= — / drf 6V{t,x,ff)W{t,x,rf) 


X [dv , ( 22 ) 

R'i 


from which we conclude 
iB f 

ydTy 5V{t,x,rf) 

R'i 

X J dvdp . (23) 

R'^xR'^ 


This equation confirms the skew-hermiticity of the ma¬ 
trix representation of the pseudo-differential operator, 
which is an important property for the stability of the 
proposed algorithm as we will see in the next section. 
An example of this matrix representation is shown in ap¬ 
pendix 1^ The entries of the matrix are given by 

^ (t>l{^)vi4>i{d) , (24) 


III. NUMERICAL METHOD 

For the numerical treatment, the expansion in Eq. 
(ED, is cut at the index A, i.e. we assume all higher 
coefficients to be zero. The problem is hence shifted to 
the time-evolution of the TV-dimensional coefficient vec¬ 
tor with the initial condition 

d{to,x) = [dv W{to,x,v)${v) , (25) 

where (j) = <(' 2 , • ■ ■, '('at)- Therefore, we work with a fi¬ 

nite set of N balance equations (PDEs) in the form of Eq. 
(ED- It is important to note that, thanks to the Cauchy- 
Kowaleski theorem, see Ref. [l^, we know that the sys¬ 
tem will locally have a unique analytical solution if the 
coefficient matrix My is an analytic function. This con¬ 
dition is sufficient, since the matrices A^l are constant. 
In addition, we would like to mention that this does not 
necessarily apply if My belongs to the larger group of 
smooth functions, see Levy’s argument in Ref. [l 8 |. 


A. Operator-splitting 


To proceed with the problem we use an operator¬ 
splitting technique (“divide-and-conquer”), i.e. we sepa¬ 
rate the action of the “streaming”, 

d 

Sd= - , (26) 

i=l 

and “forcing”, 

IFtd =—My{t,x)d , (27) 


operators to apply them sequentially. First, we discretize 
the time interval from zero to t in Nt periods of duration 
6t. Then we can write the approximated solution to Eq. 
( 1 ^ as 


Aft-l 


i{t,x) ~ exp 


k=0 


{k-\-l)St 


SSt j dt' Tt' do(^), 
k&t j 


Nt-l 

n 

k^O 


e^^* exp 


( {k+l)St 

J dt' Tt' do(x) 

Y kSt 


Nt-l 


JJ- e‘5‘5*e-^'=^“^*do(f) -b 0(St) , 


(28) 


fe =0 


where in the first step we have used the third-order ac¬ 
curate Fer expansion [l^; in the second step, simple op¬ 
erator splitting; and the numerical integration procedure 


(k-i-l)St 

J dt'Tt'- TkstSt + 0{6t‘^) , 

k6t 


which shows that it is hermitian. 


( 29 ) 
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in the third step. It is important to apply the operators in 
a time-ordered product series, which is indicated by the 
arrow above the product sign. The obtained method will 
be first-order accurate if it is stable and the numerical 
procedure for each operator (streaming and forcing) is at 
least second-order accurate. The total error arises since 
the matrices My are in general not commut¬ 

ing and because of the second-order accurate integration 
procedure. For a second-order accurate method we write 


Nt-l 


~ e‘®'^*exp 


k=0 


f {k+l)5t 

J dt' Tv I ao(®)> (30) 

y kst 


a*{t,x) = e x) , 


where we have used the Strang-splitting [ 2 ^ . To achieve 
the demanded accuracy we have to use a third-order 
accurate integration formula for the forcing operation, 
whereas second-order accuracy in the definition of a* is 
sufficient, because it acts only twice during the evolution. 
If the potential has an explicit time-dependence one can 
use the midpoint rule, 


J dt' Tv^Tt^k+^^stSt + 0{St^) . (31) 

kSt 


operator is a unitary rotation of the coefficient vector, 
whose matrix form can in general be calculated before 
starting the simulation. For a skew-symmetric matrix 
one could use the method described in Ref. or a 
Fade approximation 


gJ^kStSt 


5t 

-1 

5t 

1 + —Mv{k5t^ x) 


1 — —Mv{kSt, x) 


For the case of the Wigner-Poisson problem one needs 
to compute the product of matrix times vector at ev¬ 
ery time-step, which for instance can be efficiently done 
with the ’’Expokit” software package [1^ or using a pre¬ 
calculated explicit formula. 

One might be tempted to use explicit schemes, such as 
Euler or Runge-Kutta, to approximate the forcing. How¬ 
ever, these methods can become unstable for strongly 
changing potentials and poor temporal and spatial res¬ 
olution, which we will show for two examples by evalu¬ 
ating the amplification factor g in von Neumann’s sta¬ 
bility analysis (25j |. Consider a time-independent one¬ 
dimensional anharmonic potential and the Euler as well 
as the fourth-order Runge-Kutta method (RK4) as ap¬ 
proximations of the forcing, whose amplification factors 
are given by 


ffEuler = |1 — 5tMv{x)\2 , (33) 


For the Wigner-Poisson problem [2l| where one needs 
to determine the self-consistent electro-static potential, 
AV = ep{t,x), at every time-step, we make use of the 
fact that the forcing operation does not change the den¬ 
sity and hence the electro-static potential. Taking Eq. 
(HID and integrating by parts we can show 

/du {dtW + e[V]W) =dtp = 0 . (32) 

S'* 


4 

5RK4 = ii+ 

i=l 


[—6tMvix)Y 


(34) 


where | ... I 2 stands for the 2- or spectral-norm, such that 
Parseval’s identity is applicable. The plots for both meth¬ 
ods are shown in Eig. [TJ One observes the big amplifica¬ 
tion factor at the domain boundary, caused by the strong 
potential variation in this area, cf. Fig. dj which may 
eventually trigger a numerical instability. 


Consequently, if the numerical procedure in this step con¬ 
serves the density up to 0{St^), it will be sufficient to re¬ 
calculate the forcing operator after each streaming, which 
coincides with our time-step definition in Eq. (I30|l . The 
questions that remain to be solved are how to compute 
approximations of the operators’ actions (“stream¬ 
ing”) and (“forcing”) such that the resulting al¬ 

gorithm is stable, computationally efficient, and of the 
desired accuracy (first- or second-order). 

B. Forcing 

As it was mentioned in section III Al the matrix My 
is skew-hermitian, which means that it belongs to the 
Lie algebra of the group of unitary matrices. Depending 
on the basis choice we might also find the subgroups of 
special unitary or special orthogonal matrices if My is 
a traceless, skew-Hermitian, complex matrix or a real, 
skew-symmetric one. Hence, the action of the forcing 


C. Streaming 

The streaming can in general be achieved by meth¬ 
ods handling (non-)linear hyperbolic systems of conser¬ 
vation laws, often used in computational fluid dynam¬ 
ics, such as finite difference, volume, elements or lat¬ 
tice Boltzmann [^ . Using, for example, in d = 1 a 
flux vector splitting [T7I| . we first diagonalize the matrix 
^ 4 ( 1 ) = TDaT~^. Then we define the new coefficient 
vector 6(t, x) = T~^d(t, x) and the modified forcing term 
My = T~^MyT, such that the new system of partial 
differential equations can be written as 

dtb + DAdxb + My{t,x)b = 0 . (35) 

To simulate the action of the streaming operator, one 
can now employ the first-order accurate upwind or the 
second-order accurate Lax-Wendroff scheme [l^, since 
only has real eigenvalues, i.e. Da is a real diagonal 
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FIG. 1. Amplification factor for an anharmonic potential 
{K = 0.5) using Enler {l/5x = 100) and RK4 [l/Sx = 50) 
methods with N = 10. 



-4-3-2-101234 

X 

FIG. 2. Dimensionless anharmonic potential for K = Q (solid) 
and K = 0.5 (dashed). 


matrix. The drawbacks of the explicit methods are that 
the Courant-Friedrichs-Levy condition 


5x 


< 1 


(36) 


needs to be fulfilled for a stable simulation {conditional 
stability) and that they introduce a considerable amount 
of dissipation especially if structures with large gradients 
are streamed [l3l| . If we are employing Hermite functions 
we could - in the spirit of the lattice Boltzmann method 
[28l - l^ - use an “exact” streaming operation, which will 
mitigate the dissipative effects. For this we perform a 
discrete Hermite transform in F-space 3l| from the coef¬ 
ficient vector to the Wigner function, stream, and trans¬ 


form back to the coefficient vector. However, a more 
detailed analysis and implementation will be the topic of 
a subsequent article. 


D. Stability 


The proposed method for evolving the Wigner func¬ 
tion will be stable if the operations streaming and forc¬ 
ing are both stable. Since the action of the forcing can 
be described as a unitary rotation one should make sure 
that the numerical technique conserves this property and 
hence has an amplification factor of unity. In that re¬ 
spect, the skew-hermiticity of My, the matrix represen¬ 
tation of 0[H], is a crucial property for the stability of 
such algorithms, which we will demonstrate in appendix 
m for an asymmetric Hermite basis. This may require a 
very accurate result for the rotation matrix or the usage 
of Clifford algebras [s^ to perform the rotation. How¬ 
ever, the least computationally expensive operation is the 
direct matrix-vector multiplication 0{N‘^), whereas the 
use of the algebra will need slightly more operations (al¬ 
though the scaling is the same). For the streaming, one 
can use any stable method that handles linear advection 
equations, such as flux vector splitting [27|, Godunov, 
finite volume or Hnite element [l3|. The resulting time- 
evolution of the Wigner function will hence be stable. 


IV. SIMULATION 

For the validation of our numerical procedure we sim¬ 
ulate the time-evolution of an (an-)harmonic oscillator. 
The advantages of these examples are that, on one hand, 
we can compare with the analytical Wigner function of 
a harmonic oscillator, which is calculated as described in 
Ref. [l^. On the other hand, we can observe the effects 
of quantum corrections to the classical dynamics for an 
anharmonic potential f7anh(^ = \'mu}\<f\^ + 

0. In the case of the double well potential f7mh(^ = 
cmu}\q\^ + ™ we can observe the tunneling phe¬ 

nomenon in the Wigner formalism, since for certain pa¬ 
rameter ranges c < 0 and K > 0 the system has states 
with eigenenergies below 0 which would not allow classi¬ 
cal particles to travel from one potential minimum to the 
other. As we have described in the introduction to the 
Wigner formalism, see section |H1 we use a dimensionless 
form of the Schrodinger and Wigner equations. For our 

examples, we take I = \ —, T = — and U = hio as 
length, time, and potential scales, respectively, to find 
e = 1 and B = 1. The dimensionless time-dependent 
Schrodinger equation reads 

idt'^ = + c|Fp + K\x\^ d;* , (37) 

such that the eigenfunctions and -values of the dimen¬ 
sionless Hamilton operator at A = 0 and c = 1/2 are 
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given by 

/Com (38) 

V7r‘^/22l”ln! 

Zn = \n\ + d/2 , (39) 

where n = (ni,..., n^) is a multi-index and 

ijh»l(x') = (-l)l”lel“l', (40) 


the d-dimensional Hermite polynomial, according to Ref. 
[ 3 ^ . The dimensionless Wigner equation in differential 
form becomes 


dtW + v-V^W -2{c+2K\x\^)x-VyW + ec[K]W = 0 , 

(41) 

where 

e,[K]W = f E (42) 

|A|=3 

is the quantum correction to the “classical” dynamics of 
the particle. 


A. Basis of Hermite functions 


In our simulation, we choose Hermite functions as or¬ 
thonormal basis set in momentum-space, i.e. 




= -|«l "/2 




The proof is given in Ref. [l^ . Using the Wigner trans¬ 
form, defined by Eq. m , we can write 

a|f'(to,^) = fdv W{to,x,v)4i^^\v) , 

R‘^ 

= Jdy d>* (^to, x + ^'^d> (^to,x-^'^ (y) 

R‘^ 

gdj|fc| 

X - T ■ 

(2^)i 

One can see that the obtained coefficients are real due 
to the symmetry properties of the Hermite functions 
^^k^i~y) = (y)- In addition, we can simplify 

Eq. (1^51) by using the same property to obtain 

i^v)ki = —f dy SV{t,x,fJ)(j)^i^^{f])(j)f^{f]) . 

e jR‘i 

(46) 

Looking at this result, one can realize that My is a real, 
skew-symmetric matrix. The matrices A^') are real, sym¬ 
metric and sparse for this basis choice. They are sparse, 
because, regardless how the basis functions are ordered at 
most two entries per row or column are non-zero due to 
the recursion relation of the d-dimensional Hermite poly¬ 
nomials [s^. For further explanations on the conserva¬ 
tion and convergence properties for this basis choice, we 
refer the reader to Ref. [^, where the authors treat the 
Vlasov equation, which can be considered as the classical 
limit, e —>■ 0, of the Wigner equation. 


B. Harmonic oscillator 


where fc is a multi-index of dimension d. Hence, Eq. m 
changes to 

N 

W{t,x,v) = E ■ (44) 

|fe|=o 

The number of basis functions that is needed to simulate 
the evolution of a given state will in general depend on 
how wide the spread of the corresponding Wigner func¬ 
tion is in momentum space. However, by scaling the Her¬ 
mite functions, cf. Ref. [s^, one can significantly reduce 
N to simulate eigenstates with higher energy. In order 
to find the necessary number of basis functions for a cho¬ 
sen accuracy one needs to take a look at the variation of 
the resulting Wigner function with respect to changes in 
N. In order to find the initial coefficients, d{tQ,x), we 
use the property of the Hermite polynomials or Hermite 
functions, defined by Eq. (H51) . that they diagonalize the 
Fourier transform operator, 

f dv ' (u) = (\/^)'^jl'"l(()|f'(y) . (45) 

R‘i 


We run a simulation of an one-dimensional harmonic 
oscillator with the second order accurate method, using a 
Lax-Wendroff scheme for the streaming, a spatial resolu¬ 
tion of 6x = 1/50, and periodic boundaries at a: = ±3.5. 
The resulting matrices for the reaction-advection system 
can be calculated using formula m in appendix El As 
initial state 'L we choose a superposition between ground 
and first excited state , since the Wigner func¬ 

tion of a single eigenstate is time-independent. Thus, we 
can observe the evolution for a system whose probability 
density changes in time. In Fig. [31 we show a compar¬ 
ison between the analytical spatial probability density 
p^{t,x) = |4'(t,x)p and the probability density calcu¬ 
lated from the Wigner function with 

pw{t,x)= dv W{t,x,v) = '^ak{t,x) dv 4>k{v) . 

R ^=0 R 

The comparison shows very good agreement. However, 
the actual results for the Wigner function are more in¬ 
sightful, since they contain additional information. They 
are shown in Figs. |3]|S] together with the contour lines 
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at W = 0 and W = ±0.025. One observes a “rigid” ro¬ 
tation of the Wigner function in phase space, which is 
typical for the harmonic oscillator. This can be seen by 
using the method of characteristics for solving Eq. 6ID 
at if = 0 and c = 0.5 which leads to solving the Hamilton 
equations 

X = V , V = —X . (47) 

The period for one revolution is T = 2'k{£i — Eq) = 27r, 
which is also confirmed by the simulation in terms of the 
temporal error margin. The contour line at IT = 0 close 
to the boundaries shows patterns which are not present in 
the analytical solution. They are caused by the numerical 
error fluctuations, see Fig. El since the magnitude of the 
Wigner function in that region becomes comparable to 
the numerical error. 



-3-2-10 1 2 3 


X 

FIG. 3. Temporal evolution of the probability density for the 
harmonic potential (c = 0.5, K — 0), (solid lines) and pw 
(points) for 4' = (4/0 -|- 4'i)/\/2 using N = 16 Hermite basis 
functions. 


C. Convergence 

Based on the work in Ref. [l^, one can calculate the 
exact Wigner transform Wex of any wave function 4/ (t, x) 
expanded in Hermite functions. We will use this for¬ 
mula to calculate the Wigner transform for an eigenstate 
4/„(t,a:) of the harmonic oscillator and compare our re¬ 
sults for different numbers of basis functions N and spa¬ 
tial resolutions l/6x. It is important to note that the ex¬ 
act Wigner function of an eigenstate for RT = 0, c = 0.5 
and c? = 1 is given by Laguerre polynomials through 

Wn{x,v) = ^^L42{x^+v^)]e-^ , (48) 

TT 



FIG. 4. Wigner function for the harmonic potential (c = 0.5, 
A' = 0) at t = 0 for the superposition T = (4^0 ± 4'i)/v/2 
using N — 16 Hermite basis functions; contour lines at IT = 0 
(white) and IT = ±0.025 (black/gray). 



FIG. 5. Wigner function for the harmonic potential (c = 0.5, 
K — 0) at t — 1.25 for the superposition 4^ = (4'o ± 4'i)/%/2 
using N = 16 Hermite basis functions; contour lines at IT = 0 
(white) and W = ±0.025 (black/gray). 


which does not give a finite expansion into Hermite func¬ 
tions in V. The deviation of our results from the ana¬ 
lytical solution after one period is shown in Fig. El We 
observe that the error is of the order of 10“'* and its mag¬ 
nitude is rather homogeneously distributed. In Fig. [TU] 
we show the convergence of the second order accurate 

















FIG. 6. Wigner function for the harmonic potential (c = 0.5, 
IF = 0) at t = 2.51 for the superposition ^ = (5'o + \['i)/\/2 
using = 16 Hermite basis functions; contour lines at IF = 0 
(white) and W — ±0.025 (black/gray). 



FIG. 8. Wigner function for the harmonic potential (c = 0.5, 
A' = 0) at t = 5.01 for the superposition 3' = (3^0 ± 3'i)/v/2 
using = 16 Hermite basis functions; contour lines at IF = 0 
(white) and IF = ±0.025 (black/gray). 


since FF(i, ±5, v) ^ 0(10“^°). Looking at Fig. [TOl we ob¬ 
serve that the second order convergence can only be veri¬ 
fied for sufficiently many basis functions (here: N = 32). 
This behavior is caused by a total error that is composed 
by the discretization of time and real-space as well as the 
approximation of the Wigner function with a finite num¬ 
ber of basis functions in momentum-space. Therefore, 
we expect A to saturate for a fixed resolution and an 
increasing number of basis functions, or in the opposite 
scenario, which can be deduced from Fig. [T0]for N = 16. 


FIG. 7. Wigner function for the harmonic potential (c = 0.5, 
AT = 0) at t = 3.76 for the superposition 3* = (3*0 ± 3'i)/\/2 
using A = 16 Hermite basis functions; contour lines at IF = 0 
(white) and IF = ±0.025 (black/gray). 


method by looking at the error 

(50) 

V 

AW(Xi,Vj,t) = W(xi,Vj,t) - Wex(Xi,Vj,t) , (51) 

for periodic boundary conditions in real-space and a do¬ 
main size X G [—5,5]. The error is evaluated by choos¬ 
ing the same momentum- and space-grid. The domain 
size is chosen such that boundary effects do not signif¬ 
icantly influence the error in the convergence analysis. 



2 • 10 "'‘ 

1 ■ 10 "'‘ 
0 - 10 ° 

-1 • 10 "'‘ 
-2 • 10 "'‘ 
-3- 10"'‘ 


FIG. 9. Error of the Wigner function for the harmonic po¬ 
tential (c = 0.5, K = 0.5) at t = 6.28 for the superposition 
3- = (3-0 ± 3'i)/y2 with N = 16, <5® = 1/50. 
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FIG. 10. Convergence analysis of the harmonic Wigner func¬ 
tion after t = 2 tt with respect to Sx and Nb- 


D. Anharmonic oscillator 

For an anharmonic potential c = 0.5 and K > 0 we 
approximate the eigenstates by a superposition of 

Nf, harmonic eigenstates, i.e. 

Nb 

. (52) 

k^O 

Then we determine the coefficient vector by diag¬ 
onalizing the matrix representation of the anharmonic 
Hamilton operator. This works very well for moderate 
anharmonicities, but becomes very costly for K > 10“^, 
as can be seen in Figs. [TT] and [T^l In addition, one also 
observes that, as expected, the ground state converges 
faster than the first excited state. 

The simulation is run with the second order accurate 
method for periodic boundary conditions at a: = ±3.5 
with a spatial resolution of l/5x = 50. The result for the 
spatial probability evolution is shown in Fig. 1131 where 
we observe a good agreement with the wave function dy¬ 
namics. In Figs. fTUrrSl we show the Wigner function 
evolution together with the contour lines at IF = 0 and 
W = ±0.025. They depict a “rotation” with a smaller 
period Tan = 27r/(£’^“”^ — < 2tt. In this case it 

is not a rigid rotation, since the Wigner function gets 
compressed in position- and broadened in momentum- 
space due to the larger potential and the particle num¬ 
ber conservation. The contour line at FF = 0 close to the 
boundaries shows again the numerical error fluctuations, 
since in this region the magnitude of the Wigner function 
becomes comparable to the error, which is 0(10“^). It 
was estimated by comparing the initial and final Wigner 
function of one revolution. 



FIG. 11. Gonvergence of eigenstate coefficient vector for 
ground and first excited states of the anharmonic potential 
(c = 0.5, K — 0.001) up to double precision. 



FIG. 12. Gonvergence of eigenstate coefficient vector for 
ground and first excited states of the anharmonic potential 
(c = 0.5, K = 0.5) up to double precision. 


E. Double well potential 

To study tunneling effects, we simulate the time- 
evolution of the Wigner function for “bounded” states 
of a one-dimensional double well potential V (x) = cx"^ ± 
Kx^, cf. Fig. [HI with the second order accurate method, 
periodic boundary conditions at a; = ±4 and spatial reso¬ 
lution l/5x = 50. We call a state bounded if its eigenen- 
ergy is smaller than zero and hence below the potential 
barrier around x = 0. Taking again the Hermite basis 
to approximate the ground and first excited state as in 
Eq. (15^ . we show in Fig. [HI the evolution of the prob- 
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FIG. 13. Temporal evolution of the probability density for 
the anharmonic potential (c = 0.5, K = 0.5), (solid lines) 
and pw (points) for T = (Tg“"' + using Nh = 150, 

and = 16 Hermite basis functions. 


FIG. 15. Wigner function for the anharmonic potential 
(c = 0.5, K = 0.5) at t = 0.77 for the superposition 
^ using Nb = 150, and iV = 16 Her¬ 

mite basis functions; contour lines at IF = 0 (white) and 
IF = ±0.025 (black/gray). 



FIG. 14. Wigner function for the anharmonic potential 
(c = 0.5, K = 0.5) at t = 0 for the superposition = 
(^(“"•) _|_ using Nh = 150, and = 16 Hermite ba¬ 

sis functions; contour lines at IF = 0 (white) and IF = ±0.025 
(black/gray). 


ability density in comparison to the evolution according 
to Schrodinger’s equation. The agreement is very good. 
In Figs. [2T]|2^ one can see the evolution of the Wigner 
function for the tunneling of the state through the po¬ 
tential barrier together with contour lines at IF = 0 
and IF = ±0.025. The error during the revolution is 
at most 0(10“'^), as can be seen in Fig. [26l In ad¬ 
dition, one observes by looking at the contour line for 
IF = —0.025 the appearance of ripples and valleys in the 
front and the back of the positive quasi-probability den- 


FIG. 16. Wigner function for the anharmonic potential 
(c = 0.5, K = 0.5) at t = 1.54 for the superposition 
3- = (![<[,“") -H 4'^“"')/y2 using Nh = 150, and A^ = 16 Her¬ 
mite basis functions; contour lines at IF = 0 (white) and 
IF = ±0.025 (black/gray). 


sity during the tunneling process of the particle through 
the potential barrier, which indicate the non-classical 
behavior in the corresponding coordinate space. The 
contour line at IF = 0 close to the boundaries shows 
again the numerical error fluctuations in regions where 
the magnitude of the Wigner function becomes compa¬ 
rable to the numerical error. The period of the revolution 
Tmh = 27r/(£^™^^ — ^ 27r is much larger than the 
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by measuring the expectation values 



FIG. 17. Wigner function for the anharmonic potential 
(c = 0.5, K = 0.5) at t = 2.23 for the superposition 
^ using Nb = 150, and iV = 16 Her- 

mite basis functions; contour lines at VF = 0 (white) and 
W = ±0.025 (black/gray). 



0.3 

0.2 

0.1 

0 

- 0.1 

- 0.2 


FIG. 18. Wigner function for the anharmonic potential 
(c = 0.5, K = 0.5) at t = 3.09 for the superposition 
^ using Nb = 150, and iV = 16 Her- 

mite basis functions; contour lines at IF = 0 (white) and 
IF = ±0.025 (black/gray). 


{Ax^){Av^)^{ix-{x))^){{v-{v)r) (53) 

= {{x^)w - {x)w) {{v^)w - (w)w) > (54) 

Cov(x, v) = + xx) — {x){v)^ /{AxAv) (55) 

__ {xv)w - {x)w{v)w _ 

- {x)w^/{x‘^)w - {v)w ’ 

where {f{x,v))w = fdxdv f(x,v)W. The first quantity 
measures the well-known standard deviation of a quan¬ 
tum state in coordinate and momentum space that ful¬ 
fills Heisenberg’s uncertainty principle AxAv > f. In 
Fig. [23 we show the coordinate and momentum uncer¬ 
tainty in the form of rectangles, i.e. the width, height 
and aread correspond to Aa;, Av and AxAv, respectively. 
In that way, one can see that the standard deviation in 
position measurements mainly contributes to the uncer¬ 
tainty and its temporal change. The second quantity 
Cov(a:, v) is the covariance between the coordinate and 
momentum variable in the corresponding Wigner func¬ 
tion normalized with the standard deviations, such that 
|Cov(x,'u)| < 1. The evolution of these expectaton values 
is shown in Fig. |2H1 One observes a periodic behavior 
with T = Tm?i/2 and finds the maximum uncertainty 
AxAv exactly when the peak of the spatial probability 
density tunnels through the potential barrier in the mid¬ 
dle of the double well potential. In contrast Cov{x,v) 
behaves similar to the first temporal derivative of the 
uncertainty. 



X 


FIG. 19. Double well potential (c = —0.4, K = 0.05) and 
probability density of ground and first excited states dis¬ 
placed from 0 by = —0.310 (dashed horizontal line) 

and = —0.173 (dotted horizontal line) for Nb = 86 

Hermite basis functions. 

one for a harmonic oscillator. In addition to the proba¬ 
bility density and the Wigner function, we have also an¬ 
alyzed the spread of the Wigner function in phase space 
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-4-3-2-10 1 2 3 4 

X 

FIG. 20. Temporal evolution of the probability density for 
the double well potential (c = —0.4, K = 0.05), pip (solid 
lines) and pw (points) for the superposition 4^ = + 

using Nh = 86, and N = 32 Hermite basis func¬ 
tions. 



FIG. 22. Wigner function for the double well potential 
(c = —0.4, K = 0.05) at t = 9.16 for the superposition 
T = (4'',’"'*’ + 4'^™'‘^)/V2 using Nb = 86, and N = 32 Her¬ 
mite basis functions; contour lines at IF = 0 (white) and 
W = ±0.025 (black/gray). 



FIG. 21. Wigner function for the double well potential 
(c = —0.4, K = 0.05) at t = 0 for the superposition 
4' = (4'[,’"'*^ ± 4'^’"'‘^)/V^ using Nb = 86, and iV = 32 Her¬ 
mite basis functions; contour lines at IF = 0 (white) and 
IF = ±0.025 (black/gray). 


V. CONCLUSIONS 

We have developed a semi-spectral simulation method 
for the time-evolution of the Wigner quasi-probability 
distribution that uses a spectral-decomposition of the 
distribution into arbitrary basis functions of in 

momentum-space, which transforms the original partial 
differential equation into an infinite-dimensional set of 


FIG. 23. Wigner function for the double well potential 
(c = —0.4, K = 0.05) at t = 18.3 for the superposition 
^ ± 4'^’"'‘^)/^/2 using Nb = 86, and iV = 32 Her¬ 

mite basis functions; contour lines at IF = 0 (white) and 
IF = ±0.025 (black/gray). 


advection-reaction equations. 

For the numerical treatment, we introduce a cut¬ 
off in the expansion, which makes the system finite¬ 
dimensional, and split the operators for the reaction and 
advection part so as to apply them sequentially to the 
distribution function. We demonstrated that, due to the 
skew-hermitian symmetry of the matrix representation 
of the pseudo-differential operator (Lie algebra), the ac¬ 
tion of the forcing or reaction operator (Lie group) is 
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1 ■ 10 "'‘ 
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FIG. 24. Wigner function for the double well potential 
(c = —0.4, K = 0.05) at t = 27.5 for the superposition 
4' = using Nb = 86, and iV = 32 Her- 

mite basis functions; contour lines at IF = 0 (white) and 
IF = ±0.025 (black/gray). 



FIG. 25. Wigner function for the double well potential 
(c = —0.4, K = 0.05) at t = 36.7 for the superposition 
4- = (4'[,"*'*’ ± 4'^”"'‘^)/^ using Nb = 86, and iV = 32 Her- 
mite basis functions; contour lines at IF = 0 (white) and 
IF = ±0.025 (black/gray). 


a unitary rotation, which stabilizes the simulation even 
for strongly varying potentials compared to other ex¬ 
plicit methods, such as Euler or RK4. The advection 
or streaming part can be handled by many numerical ap¬ 
proaches from computational fluid dynamics. Here, we 
have chosen a flux-vector splitting for the validation of 
our method by simulating a single, non-relativistic, spin¬ 
less particle subject to a one-dimensional (an-)harmonic 
or double well potential with Hermite basis functions. 


FIG. 26. Error of the simulated Wigner function for the dou¬ 
ble well potential (c = —0.4, K = 0.05) at t = 45.9 for the 
superposition 4^ = (4 'q™'^^ ± 4'5^™'^^)/v/2 using Nb = 86, and 
N = 32 Hermite basis functions. 



FIG. 27. Temporal evolution of the spread in the form of rect¬ 
angles (width = Ax, height = Av, area = AxAv) of the simu¬ 
lated Wigner function for the double well potential (c = —0.4, 
K = 0.05) for the superposition 4t = (4tg"‘^^ ± 4'J"‘^^)/%/2 us¬ 
ing Nb = 86, and N = 32 Hermite basis functions. 


Having the exact Wigner function of the harmonic os¬ 
cillator, we verified the second-order convergence of the 
method and also demonstrated its applicability to non- 
classical dynamics in the case of strong anharmonicities 
and tunneling phenomena. 

The disadvantage of an arbitary basis choice is 
the higher computational cost of 0{N‘^) compared 
to 0{N log N) for a Fourier basis, since the pseudo¬ 
differential operator is diagonal for this basis choice, as 
shown in Refs. 0, 0 . However, if one only considers 
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FIG. 28. Temporal evolution of the uncertainty (solid) and 
covariance (dashed) of the simulated Wigner function for the 
double well potential (c = —0.4, K = 0.05) for the superpo¬ 
sition ^ -h using Nb = 86, and N ^ 32 

Hermite basis functions. 


momentum-derivatives up to order Nx <C N and an ex¬ 
plicit scheme such as fourth order Runge-Kutta is used, 
the computational cost also scales like 0{N) [s^. In 
addition, the artificial periodization of the Wigner dis¬ 


tribution in momentum-space, caused by the plane wave 
approximation, lives in a different function space than 
the original Wigner function, thus giving rise to unphysi¬ 
cal self-interactions at the domain boundaries [l^ . These 
basis functions are also not well suited to the simulation 
of structures that are strongly localized in momentum- 
space, such as particles in periodic potentials, since this 
would require a very large number of such functions. 

The CPU time for the time evolution of one revolu¬ 
tion for the harmonic Wigner function, using the second- 
order accurate method with N = 32 Hermite basis func¬ 
tions and a spatial resolution l/5x = 100, i.e. 700 grid 
points and 4558 time-steps, is approximately 7.92 s using 
a single core of a 3 GHz Intel(R) Core(TM)2 Quad CPU 
Q9650 processor. 

As future work, we plan to study phase transitions in 
open quantum systems, the effects of scattering quan¬ 
tum Boltzmann equation"), for example in the case of 
electrons and phonons in semiconductor devices, the ef¬ 
fects of boundary conditions, cf. Ref. and stochastic 
perturbations. Furthermore, we plan to analyse the in¬ 
fluence of decoherence on the topology of the Wigner 
function in two dimensions, cf. Ref. |35l |. 
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Appendix A: Matrix-representation of 
pseudo-differential operator 


For a one-dimensional, analytical potential and the 
Hermite function basis the matrix-representation of the 
pseudo-differential operator simplifies from Eq. (1461) to 

OO o 

_ / f \ 2n 

= Mr,dl-+^V{t,x) , (Al) 

n—0 

p Nin+l 

(M„)fc./ = J d?7 ■ 

R 


Looking at this result one observes how the contributions 
from odd higher order derivatives scale with the effective 
Planck constant and the change in sign. The examples 
for A = 5 shows the filling of higher order matrices with 
more and more entries. In the case of A = 5, the matrix 
Mv is already filled for n = 2, i.e. considering the fifth 
derivative of the potential. As soon as My is completely 
filled an explicit method, such as Euler or Runge-Kutta 
will have a computational complexity of 0{N‘^), although 
the prefactor will be smaller than for the method which 
uses the corresponding rotation matrix. 
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Appendix B: Example for unstable “basis” choice 

Assuming we are dealing with a harmonic potential, 
then the quantum corrections vanish and the Wigner and 
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Vlasov equation are identical. Using an asymmetric Her- 
mite basis, as described in [Sa, i.e. 


„2 N 


= (Bl) 


the matrix representation of 0[U] will be lower triangu¬ 
lar, which can be seen using formula m and integration 
by parts. For = 4 we find 


Mvix) = 


/ 0 0 0 0 0 \ 

0 0 0 0 
0 2 0 0 0 

0 0 V6 0 0 

V 0 0 0 2^2 0 


which is not skew-hermitian or -symmetric anymore. Ex¬ 
amining the resulting exact forcing action, we find 


^-Mv{x)St _ 


( 1 
\/2x^(5t^ 
_U3 


0 

1 

2x6t 

\/6x^5t^ 

4x^St^ 

^/3 


0 
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'/hxSt 

2\/3x2(5t^ 


0 0 \ 

0 0 

0 0 

1 0 

2^f2xbi ij 


from which we conclude for the amplification factor 


5AS = 


Mv {x)5t 


>lif(5t>0,a;^0 


This means that the method will become unstable at a 
certain time in the evolution, as described in Ref. [s^. 
There are ways to tackle t his p roblem by introducing a 
collision operator, see Re f. 13711 . or a filtering technique, 
further explained in Ref. [38|. 















